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Abstract 

Dynamical systems that describe the escape from the basins of attraction of stable in- 
variant sets are presented and analyzed. It is shown that the stable fixed points of such 
dynamical systems are the index-1 saddle points. Generalizations to high index saddle 
points are discussed. Both gradient and non-gradient systems are considered. Preliminary 
results on the nature of the dynamical behavior are presented. 

1 The gentlest ascent dynamics 

Given an energy function V on R", the simplest form of the steepest decent dynamics (SDD) 
associated with V is 

x = -W(x). (1) 

It is easy to see that if x(-) is a solution to ([T|), then U(x(f)) is a decreasing function of t. 
Furthermore, the stable fixed points of the dynamics ([T]) are the local minima of V. Each local 
minimum has an associated basin of attraction which consists of all the initial conditions from 
which the dynamics described by ((T|) converges to that local minimum as time goes to infinity. 
For ([1]), these are simply the potential wells of V. The basins of attraction are separated by 
separatrices, on which the dynamics converges to saddle points. 

We are interested in the opposite dynamics: The dynamics of escaping a basin of attraction. 
The most naive suggestion is to just reverse the sign in ([1}, the dynamics would then find the 
local maxima of V instead. This is not what we are interested in. We are interested in the 
gentlest way in which the dynamics climb out of the basin of attraction. Intuitively, it is clear 
that what we need is a dynamics that converges to the indcx-1 saddle points of V. Such a 
problem is of general interest to the study of noise-induced transition between metastable states 
[31 [5] : Under the influence of small noise, with high probability, the escape pathway has to go 
through the neighborhood of a saddle point [5]. 
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The following dynamics serves the purpose: 



x = - VV(x) + (2a) 
v = -VV(x)v + ^^v. (2b) 

We will show later that the stable fixed points of this dynamics are precisely the index- 1 
saddle points of V and the unstable directions of V at the saddle points. Intuitively the idea is 
quite simple: The second equation in attempts to find the direction that corresponds to the 
smallest eigenvalue of V 2 V% and the last term in the first equation makes this direction an ascent 
direction. 

This consideration is not limited to the so-called "gradient systems" such as ([1]). It can be 
extended to non-gradient systems. Consider the following dynamical system: 

x = F(x). (3) 

We can also speak about the stable invariant sets of this system, and escaping basins of attraction 
of the stable invariant sets. In particular, we can also think about finding index-1 saddle points, 
though in this case, there is no guarantee that under the influence of small noise, escaping the 
basin of attraction has to proceed via saddle points [9]. 
For non-gradient systems, ([2]) has to be modified to 



(w,v) 

v = (VF(x))v - a(v)v, 

w = (VF(x)) T w - /3(v, w)w. 



(4a) 

(4b) 
(4c) 



Here two directional vectors v and w are needed in order to follow both the right and left 
eigenvectors of the Jacobian. Given the matrix VF(x), two scalar valued functions a and j3 are 
defined by 
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a(v) = (v,(VF(x))v), 



(5a) 



/3(v, w) = 2(w, (VF(x))v) - a(v). 



(5b) 



We have taken and we will take the normalization such that (v, v) = 1 and (w, v) = 1. They are 
to keep the normalization such that (v, v) = 1 and (w,v) = 1. This normalization is preserved 
by the dynamics as long as it holds initially. Thus, the first equation in (j4|) actually is equivalent 
to x = F(x) — 2(F(x),w)v. (Of course, one can enforce other types of normalization condition, 
such as the symmetric one: (v, v) = (w, w) and (w, v) = 1, and define new expressions of a and 
P accordingly.) In the case of gradient flows, we can take w = v and ^ reduces to ©. 



Figure 1: Illustration of the gentlest ascent dynamics. F is the force of the original dynamics 
and F is the force of the gentlest ascent dynamics, vi and V2 represent the unstable and stable 
right eigenvectors, respectively; Wi and W2 are the corresponding left eigenvectors. Note that 
wi _L V2 and W2 _L vi. F has the decomposition F = Fi +F2 = civi +C2V2 where the coefficient 
ci = (F, wi)/(vi,wi). Thus, F := -Fi + F 2 = F - 2Fi = F - 2ciVi. 

We call this the gentlest ascent dynamics, abbreviated GAD. It has its origin in some of 
the numerical techniques proposed for finding saddle points. For example, there is indeed a 
numerical algorithm proposed by Crippen and Scheraga called the "gentlest ascent method" [2] ■ 
The main idea is similar to that of GAD, namely to find the right direction, the direction of the 
eigenvector corresponding to the smallest eigenvalue and making that an ascent direction. But 
the details of the gentlest ascent method seem to be quite a bit more complex. The "eigenvector 
following method" proposed in literature, for example, [TJ [3], is based on a very similar idea. 
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There at each step, one finds the eigenvectors of the Hessian matrix of the potential. Also closely 
related is the "dimer method" in which two states connected by a small line segment are evolved 
simultaneously in order to find the saddle point [TJ. One advantage of the dimer method is that 
it avoids computing the Hessian of the potential. From the viewpoint of our GAD, the spirit 
of "dimer method" is equivalent to use central difference scheme to numerically calculate the 
matrix-vector multiplication in GAD ((4]) and ([5]) by writing (VF(x))b = ^F(x + eb)| e= o ~ 
^(F(x + eb) — F(x — eb)) for any vector b. 

We believe that as a dynamical system, the continuous formulation embodied in ([2]) and ((4]) 
has its own interest. We will demonstrate some of these interesting aspects in this note. 
Proposition. Assume that the vector field F is C 3 (R"). 

(a) If (x,,v,,w») is a fixed point of the gentlest ascent dynamics Q and v*, w* are normalized 
such that v!f v* = v^w» = 1, then v» and w* are the right and left eigenvectors , respectively, 
of VF(xh.) corresponding to one eigenvalue A*, i.e., 

(VF(x,))v, = A*v», (VF(x,)) T w, = A*w*, 

and x* is a fixed point of the original dynamics system, i.e., F(x*) = 0. 

(b) Let x s be a fixed point of the original dynamical system x = F(x). If the Jacobian matrix 
J(x s ) = VF(x s ) has n distinct real eigenvalues Ai,A2,--- , A„ and n linearly independent 
right and left eigenvectors, denoted by v,; and Wj correspondingly, i.e., 

J(x s )vi = AjVj, JT(x s ) T Wi = AjWi, i = l,---,n 

and in addition, we impose the normalization condition v^v,; = w^v.; = 1, Vi , then for all 
i = 1, ■ ■ ■ , n, (x s , Vj, Wj) is a fixed point of the gentlest ascent dynamics (jj). Furthermore, 
among these n fixed points, there exists one fixed point (x s ,Vj/,Wj') which is linearly stable 
if and only if x s is an index- 1 saddle point of the original dynamical system x = F(x) and 
the eigenvalue \i> corresponding to v^, Wj< is the only positive eigenvalue of J(x s ). 

Proof. 

(a) Under the given condition, it is obvious that (VF(x,))v, = a(v*)v» and (VF(x t )) T w t = 
/3(v*, w»)w*. By definition and other conditions, /3(v*,w„) = 2w^(VF(x»))v, — a(v„) = 
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2w^(a(xH<))v* — a(v*) = a(v*). Therefore, and w* share the same eigenvalue A* = a(v*) = 
/3(v*, w*). From the fixed point condition F(x») — 2(w;f F(x*))v* = 0, we take the inner product 
of this equation with w* to get w^F(x„) = 2w^F(x„). So w!fF(x*) = and in consequence, 
the conclusion F(x») = holds from the fixed point condition F(x») — 2(w;f F(x»))v» — again. 

(b) It is obvious that for all i, (x s , Vj, w,;) is a fixed point of the gentlest ascent dynamics ((4]) 
by the definition of Vj and Wj. It is going to be shown that we can explicitly write down the 
eigenvalues and eigenvectors of GAD at any fixed point (x s ,Vi,Wj). 

Let J(x) = VF(x). The Jacobian matrix of the gentlest ascent dynamics ((4]) has the following 
expression: J(x, v, w) = 

' (I - 2vw T )J(x), -2(F(x),w)I, -2vF(x) T 

L x , J(x)-a(v)I-vv T (J(x)+J(x) T ), 

v L 2 , -2ww T J(x)+wv T (JJ(x)+J(x) T ), JJ(x) T -/3(v,w)I-2wv T J(x) T y 

(6) 

where Li , L 2 arc n x n matrices and I is the n x n identity matrix. To derive the above formula, 
we have used the results from © that V v (a) = v T (J T + J), V v (/3) = 2w T J - v T (JJ T + JJ) and 
VwCS) = 2v T J T . 

In the first n rows of J, there are two nxn blocks which contain the term F(x) and thus vanish 
at the fixed point x s . So the eigenvalues of J(x s , Vj, Wj) can be obtained from the eigenvalues of 
its three nxn diagonal blocks: N, M and K: 

N = (I-2v,wf)J(x s ), 

M = J(x s ) - Ail - v,vf (J(x s ) + A^), 

K = J T (x s ) - Ail - 2w. t vf J T (x s ). 

Here the obvious facts that a(vi) = /3(vi,Wi) = Xi and vj I T = XivJ are applied. 

Now we derive the eigenvalues of N, M and IK by constructing the corresponding eigenvectors. 
Note that vfwj = Sij holds under our assumption of the eigenvectors. One can verify that 

Nv, = (I-2viwf)AiVi = -A,v 4 , 
Mvi = -2AiVivf v, = -2AiV^, 
Iw, = -2AjWjvf w, = -2AjWi, 
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and for all j =/= i 



(I-2v i wf)A j v j = A j v J - ! 



(7) 



(Xj — Xi)wj — 2Xj\v i v i w 



(Xj - X l )w j , 



(8) 



and with a bit more effort, 



M ( v j - (vfvjOvi) = M Vj - vf Vj(Mvi) = Mvj + 2A*(vf v 3 -)v< 

= (Aj - Ai)vi - (Aj + Ai) Vi (vf Vj ) + 2A, (vf v,)v, 
= (Aj - Ai)(v 3 - - (vfvj-)vi). 

Hence the eigenvalues of the Jacobian J at any fixed point (x s , Vj, w,) (i = 1, ■ • ■ , n) are 



The first and last set of eigenvalues have multiplicity 2. The linear stability condition is that all 
numbers in (jf 0[) are negative. Thus one fixed point (x s , vv, wv) is linearly stable if and only if 
Xi> > and all other eigenvalues Xj < for j ^ i' , in which case the fixed point x s is index- f 
saddle. 

Next, we discuss some examples of GAD. 

Consider first the case of a gradient system with V(x) = x T Ax/ (x T x), where A is a symmetric 
matrix. V is nothing but the Raylcigh quotient. A simple computation shows that the GAD for 
this system is given by: 



I v = - Aw + (v T iv)v. 

Next, we consider an infinite dimensional example. The potential energy functional is the 
Ginzburg-Landau energy for scalar fields: I(u) = J Q (i|Vu| 2 + \(u 2 — f) 2 ) dx. The steepest 
decent dynamics in this case is described by the well-known Allcn-Cahn equation: 



2A l; -Aj, {Aj : j ^ i}, {Xj - X, : j ^ i}. 



(10) 




(11) 



dtu = Au — (u 2 — l)u. 



(12) 
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A direct calculation gives the GAD in this case: 



dtu = Au — (u 2 



l)u-2(Au - (u 2 



l)u, v)v, 



< 



(13) 



d t v = Av — (3u 2 



l)v - {Av - (3u 2 



l)v,v)v, 



where the inner product is defined to be: 



(u,v) 




u(x)v(x)dx. 



Clearly both the SDD and the GAD depend on the choice of the metric, the inner product. 
If we use instead the _ff _1 metric, then the SDD becomes the Cahn-Hilliard equation and the 
GAD changes accordingly. 

2 High index saddle points 

GAD can also be extended to the case of finding high index saddle points. We will discuss how 
to generalize it to index-2 saddle points here. There are two possibilities: Either the Jacobian JJ 
at the saddle point has one pair of conjugate complex eigenvalues or it has two real eigenvalues 
at the saddle point. We discuss each separately. 

Intuitively, the picture is as follows. We need to find the projection of the flow, F(x), on the 
tangent plane, say P, of the two dimensional unstable manifold of the saddle point, and change 
the direction of the flow on that tangent plane. For this purpose, we need to find the vectors vi 
and V2 that span P. In the first case, we assume that the unstable eigenvalues at the saddle point 
are Ai,2 = A# ; ± zAj. In this case there are no real eigenvectors corresponding to Ai2- However, 
for any vector v in P, (VF)v simply rotates v inside P. Hence, V2 can be taken as (VF)vi if we 
have already found some vi £ P. The latter can be accomplished using the original dynamics in 



Using the fact that the eigen-plane of (VF) T corresponding to Aij ± iXi, which is spanned by wi 



To see how one should modify the flow F on the tangent plane, we write 




i>2 
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and w"2 = (VF) T Wi, is orthogonal to Vj for all j > 2, we can derive a linear system for c\ and 
C2 by taking the inner product of F and Wi ,w 2 . The solution of that linear system is given by: 



&22/1 — ^12/2 an/2 — 021/1 

ci = , c 2 



0-11022 — 021021 



O11O22 — a2lfl21 



(14) 



where = (w^v-,) and fj = (F(x),Wj) for i,j = 1,2. The gentlest ascent dynamics for the x 
component is 

F = F 2civi - 2c 2 v 2 . 
To summarize, we obtain the following dynamical system: 



(15) 



x = F - 2ciV! - 2c 2 v 2 , 
V! = (VF(x))vi - a(vi)vi, 
wi = (VF(x)) T wi - /3(vi,wi)wi, 
v 2 = VF(x)vi, 
_ w 2 = (VF(x)) T w l5 

where ci,C2 are given by (fT4")l and a, (3 are defined by ©• 

If the Jacobian has two positive real eigenvalues at the saddle point, say, Ai > X2 > > A3 > 
• • • , let us define a new matrix by the method of deflation: 



h 



VF 



(vi,(VF)vi) T 
7 u r v i w i ■ 

(Vi, Vi)(Wi, Vi) 



(16) 



It is not difficult to see that if Vi is an eigenvector of VF corresponding to Ai, then J2 shares the 
same eigenvectors as J, and the eigenvalues of J2 become 0, A2, A3, • • ■ . The largest eigenvalue of 
J 2 at the index-2 saddle point becomes A 2 . One can then use the dynamics (|4b| associated with 
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the new matrix J2 to find V2. Therefore, we obtain the following indcx-2 GAD 

x = F 2civi - 2c 2 v 2 , 
vi = (VF(x))vj - aivi, 
< wi = (VF(x)) T Wl -ftwi, (17) 

v 2 = JJ 2 v 2 - a 2 v 2 , 
^ w 2 = J 2 r v 2 - /3 2 w 2 , 

with the initial normalization condition (vi,vi) = (v 2 ,v 2 ) = (wi,vi) = (w 2 ,v 2 ) = 1. ci and c 2 
are given in the same way as shown above (|14[) and ai, 2 , /?i j2 are defined as follows to enforce that 
the normalization condition is preserved : ol\ = (vi, (VF(x))vi), /3i = 2(wi, (VF(x))vi) — ax 
and a 2 = (v 2 , JJ 2 v 2 ), /3 2 = 2(w 2 , J 2 v 2 ) - a 2 . 

The generalization to higher index saddle points with real eigenvalues is obvious. 

3 Examples 

3.1 Analysis of a gradient system 

To better understand the dynamics of GAD, let us consider the case when a different relaxation 
parameter is used for the direction v: 

x = -W(x) + 2(V7,v)v, 
tv = -V 2 F(x)v + (v,VVv)v. 

To simplify the discussions, we consider the limit as r — > 0. In this case, we obtain a closed 
system for x: 

x = -W(x)+2(VV,v(x))v(x), (18) 

where v(x) is the eigenvector of V 2 F(x) associated with the smallest eigenvalue. Now we consider 
the following two dimensional system: 
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where [i is a positive parameter. x± = (±1, 0) are two stable fixed points and (0, 0) is the indcx-1 
saddle point. The eigenvalues and eigenvectors of the Hessian at a point x = (x,y) are 



Ai = 3x 2 — 1 and vi = (1,0), 



A2 = /x and v 2 



(0,1). 



Therefore, the eigendirection picked by GAD is 



Vl, 



if \x\ < 




(19) 



vgad(x) 



v 2 , 




Consequently, by defining 



Vx(x t y) 




I) 2 + 2^ 2 , 



and 




i) 2 - \m\ 



we can write the gentlest ascent dynamics (|18p in the form of a gradient system driven by the 
new potential: 



where l.(x) is the indicator function. Note that Vqad is not continuous at the lines x = ± 



The point (0,0) becomes the unique local minimum of Vi, with the basin of attraction {(x,y) : 
— l<x< 1}. Outside of this basin of attraction, the flow goes to (x = ±00, y = 0) and the 
potential V\ falls to —00 . For V2, the point (0,0) is the unique local maximum and all solutions 
go to (x = ±1, y = ±00). 

If we start the gentlest ascent dynamics with the initial value x± = (±1,0), then there are 
two different situations according to whether [i > 2 or /i < 2. Although x± becomes a saddle 
point for any /z ^ 2, the unstable direction for fj, < 2 is ±v 2 while the unstable direction for 
fi > 2 is ±Vi, as illustrated in figure [3] Furthermore, from figure [3] and the above discussion, 
it is clear that the basin of attraction of the point (0, 0) associated with the potential Vqad is 
the region — ^Mp < x < \J^^- for /i < 2 and — 1 < x < 1 for fj, > 2. (which is larger than 
the basin of attraction for the Newton-Raphson method, confirmed by numerical calculation.) 




(20) 
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Figure 2: The discontinuity of VaAD(x,y = 0) at the location x = ±y Left: \i < 2; Right: 
H > 2. 



Consequently, the GAD with an initial value (xq, yo) near the local minimum x± of V converges 
to the point (0, 0) of our interest when // > 2 and |xo| < 1. 

This discuss suggests that GAD may not necessarily converge globally and instabilities can 
occur when GAD is used as a numerical algorithm. When instabilities do occur, one may simply 
reinitialize the initial position or the direction. 

3.2 Lorenz system 

Consider 

x = <r(y- x), 

V = px-y- xz, 

z = —j3z + xy. 



(21) 



The parameters we use are a = 10, /3 = | and (3 — 30. There are three fixed points: the origin 
O = (0, 0, 0) and two symmetric fixed points 



= (±y/0(p - 1), ± V£0 -!),/>- !)• 



O is an index- 1 saddle point. The Jacobian at Q± has one pair of complex conjugate eigenvalues 
with positive real part. In our calculation, we prepare the initial directions vo and Wo by running 
the GAD for long time starting from random initial conditions for v and w while keeping x fixed, 
although this is not entirely necessary. Figure |4] shows two solutions of GAD. For the index- 1 
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saddle point O, figure [5] depicts how the trajectory of GAD converges to it. It can be seen 
that the component of the original force F along the unstable direction of O is nearly projected 
out, thus the trajectory will not be affected by the unstable flow in that direction and avoids 
departing the saddle point. Therefore the trajectory tends to follow the stable manifold toward 
the saddle point when the trajectory is close enough to the saddle point. Similar behavior is 
seen for the case of searching the point Q+ which has one pair of complex eigenvalues. The 
trajectory surrounding Q + in the figure [4] spirals to Q+ an d these spirals are closer and closer 
to the unstable manifold of in the original Lorenz dynamics, which looks like a twisted disk. 
The convergence rate of the spiraling trajectories in GAD is very slow because the real part of the 
complex eigenvalues (A = 0.1474 ± 10.5243 i) in the original dynamics is rather small compared 
with its imaginary part. 




Figure 4: The trajectories of GAD for the Lorenz system starting from two initial points. They 
converge to the index-2 saddle point Q+ (marked by the dot) and the index-1 saddle point O 
(marked by "+" ) respectively. 



If we reverse time t — > —t, we have the time-reversed Lorenz system, in which the origin O 
becomes an index-2 saddle point. We can apply the index-2 GAD algorithm (|17l) to search for 
this saddle point. The GAD trajectory in this case is also plotted in the figure OH It is similar to 
the situation of GAD applied to the original Lorenz system in the sense that the GAD trajectory 
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nearly follows the z axis when approaching the limit point O. Indeed, as far as the x-componcnt 
is concerned, the linearized gentlest ascent dynamics for the original Lorenz system and the 
time-reversed one are the same. From the proof of the Proposition (particularly, note that the 
eigenvalues of N are — Xi and Xj), it is not hard to see that the eigenvalues of the linearized 
gentlest ascent dynamics at the point O are all negative and have the same absolute values as 
the eigenvalues of the original dynamics, and the two dynamics share the same eigenvectors 
(again, we mean the x component of the GAD). Thus, since the change t — > — t does not change 
the absolute values of the eigenvalues of the original dynamics, the gentlest ascent dynamics 
for the original and time reversed Lorenz system have the same eigenvalues: Ai = —23.3955, 
A2 = —2.6667, A3 = —12.3955. The two linearized GAD flows near the point O are the same: 
x(i) = e -23.3955t Vi _|_ e -2.6667t V2 + e -i2.3995t V3 ^ w h ere vi j2 ,3 are the eigenvectors: v 2 = (0, 0, 1), 

and Vi, V3 are in the z — plane. As t — s- +00, we then have x(i) ~ e _2 " 6667t V2. This explains 
why both trajectories in the figure [5] follow the z axis when approaching the saddle point O. 




0.3- I 

0.2- 

0.1 - 

-1 -0.5 0.5 1 

X 

Figure 5: How the GAD trajectories approaches the saddle point O. The curve with two arrows 
is the trajectory of index- 1 GAD for the Lorenz system; the curve with single arrow is the 
trajectory of index-2 GAD for the time reversed Lorenz system. The unstable manifold of O, 
which is tangent to the z = plane, is also shown. 
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3.3 A PDE example with nucleation 

Let us consider the following reaction-diffusion system on the domain x G [0, 1] with periodic 
boundary condition: 



The parameter <5 is fixed at 0.01 and we allow the parameter \x to vary. There are two stable 
(spatially homogeneous) solutions for certain range of fi: u + = (u+, v+) and = (0, 0). If one uses 
the square-pulse shape function as a initial guess in the Newton- Raphson method, no convergence 
can be achieved in most situations. We applied the index-1 GAD method to this example. The 
initial conditions for GAD are constructed by adding a small amount of perturbations around 
either stable solutions: u + or 0. We observed that for a fixed value of /i, the solutions of GAD 
constructed this way converge to the same saddle point. The different saddle points obtained 
from GAD at different values of \x are plotted in figure O It is also numerically confirmed that 
these saddle points indeed have index 1 and the unstable manifold goes to u + in one unstable 
direction and to in the opposite unstable direction. It is interesting to observe the dependence 
of the saddle point on the parameter fj, and that such a dependence is highly sensitive when /i is 
close to —1.046 ~ —1.045. In fact, there exists a critical value n* in this narrow interval at which 
the spatially extended system (|22[) has a subcritical bifurcation, which does not appear in the 
corresponding ODE system without spatial dependence. We refer to [I] for further discussions 
about this point. 

4 Concluding remarks 

We expect that GAD is particularly useful for handling high dimensional system in the sense 
that it should have a larger basin of attraction for finding saddle points, than, for example, the 
Newton-Raphson method. There arc many questions one can ask about GAD. One question is 
the convergence of GAD as time goes to infinity. Our preliminary result shows that GAD does 




(22) 



where 




(u — u 3 + 1.2)v + \fiu, 



(23) 
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Figure 6: The profiles of saddle points of the example (j2"3"]l (6 = 0.01). Only the component u 
is shown since v = \u 2 at the saddle point. From inside to outside, the values of fi are —1.000, 
-1.020, -1.040, -1.045, -1.046, -1.050. 



not have to converge. For finite dimensional systems, there is always local convergence near the 
saddle point. The situation for infinite dimensional systems, i.e. PDEs, seems to be much more 
subtle. Another interesting point is whether one can accelerate GAD. For the problem of finding 
local minima, many numerical algorithms have been proposed and they promise to have much 
faster convergence than SDD. It is natural to ask whether analogous ideas can also be found for 
saddle points. 
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